Quantum phase transition in space in a ferromagnetic spin-1 Bose-Einstein condensate 
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I. INTRODUCTION 



Bogdan Damski and Wojciech H. Zurek 
Theoretical Division, Los Alamos National Laboratory, MS-B213, Los Alamos, NM 87545, USA 

A quantum phase transition between the symmetric (polar) phase and the phase with broken 
symmetry can be induced in a ferromagnetic spin-1 Bose-Einstein condensate in space (rather than 
in time). We consider such a phase transition and show that the transition region in the vicinity 
'itif/ . of the critical point exhibits scalings that reflect a compromise between the rate at which the 

ISL ' transition is imposed (i.e., the gradient of the control parameter) and the scaling of the divergent 

' healing length in the critical region. Our results suggest a method for the direct measurement of 

, the scaling exponent v. 

^ '. 
Ph. 

<: 

' Studies of phase transitions have traditionally focused on equilibrium scalings of various properties near the crit- 
I ical point of a homogeneous system. The dynamics of phase transitions presents new interdisciplinary challenges. 
I— I Nonequilibrium phase transitions mayplay a role in the evolution of the early Universe [l| . Their analogues can be 
' studied in condensed matter systems Q. The latter observation led to the series of beautiful experiments 3 (see [3| 
^ , for an up-to-date review) and to the development of the theory based on the universality of critical behavior . The 

a recent progress in cold atom experiments allows for the temporal and spatial control of different systems undergoing 
I a quantum phase transition (QPT) 0, Q. These experimental developments call for an in-depth understanding of 
non-equilibrium QPTs. 

^ _ A QPT is a fundamental change in the ground state of the system as a result of small variations of an external 
. parameter, e.g., a magnetic field In contrast to thermodynamic phase transitions, it takes place ideally at 
temperature of absolute zero. 

The problem of how a quantum system undergoes a transition from one quantum phase to another due to time- 
dependent (temporal) driving has attracted lots of attention lately 8, 9j, [l^, |Tl|, [l^, [13, [3j [3 ■ Basic insights into the 
QPT dynamics can be obtained through the quantum version [8j,|l6| of the Kibble-Zurek mechanism (KZM) [3,[2]. 
(-H ' The KZM recognizes that the time evolution of the quantum system is adiabatic far away from the critical point where 
Q ' the gap in the excitation spectrum is large. The system adjusts to all changes imposed on its Hamiltonian. Near 
O i the critical point, however, the gap closes precluding the adiabatic evolution and resulting in the non-equilibrium 
dynamics. This switch of behavior happens when the system reaction time h/ A (A is excitation gap) becomes 
. comparable to the time scale on which the critical point is approached, e/\de/dt\ {e = \q — qd, q is the parameter 
^ ' driving the transition, qc is the location of the critical point). This brings us to 



(1) 



OO ■ A \de/dt\' 

, To solve ([T]) we define 1^9(^)1 = '''q^- The solution of ^ gives us time i, left to reaching the critical point, when the 
' non-equilibrium dynamics starts 
00 ' 

(ft/Ao)i/(i+"'^V^"/(i+"''\ (2) 



Above we have set A — AqI^c — (zl^" near the critical point (z and v are critical exponents). 

The quantum KZM proposes that the non-equilibrium evolution in the neighborhood of the critical point is to 
$H ' a first approximation impulse (diabatic): the state of the system does not change there. Suppose that the system 
evolves towards the critical point from the ground state far away from it. Thus its properties near the critical point 
are determined by its ground state properties at a distance tJrQ from the critical point (the adiabatic regime was 
abandoned there). As have been shown in spin models (e.g., [1,[3]) and ferromagnetic spin-1 Bose-Einstein condensates 
[ill , [l^ . [isj , this simplification allows for a correct analytical computation of the density of excitations resulting from 
the non-equilibrium quench. 

In this paper we study a quantum phase transition induced by spacial rather than temporal driving. It means that 
the driving parameter, q, depends on position and is independent of time. Such a transition was recently studied in 
the quantum Ising model [o, [ll] and the mean-field Ginzburg-Landau theory [l3| • Here we present the theory of the 
spatial quench in a ferromagnetic spin-1 Bose-Einstein condensate, which is one of the most flexible physical systems 
for studies of QPTs. 

We focus on the ground state. In the simplest approximation, the local homogeneous approximation (LHA), the 
system is locally characterized by its homogeneous ground state properties perfectly tracking spatial variations of q. 
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This approximation is good far enough from the critical point, where the heahng length [l9| is small compared to the 
imposed scale of spatial driving. Using the time quench analogy, the system evolves perfectly adiabatically in space 
away from the spatial critical point, Tc, where q{rc) — Qc- 

Around the critical point, however, the healing length diverges and so the LHA breaks down. The system enters 
"non-equilibrium" regime similarly as during a time quench. The spatial coupling term in the Hamiltonian (a gradient 
term in mean- field theories, a spin-spin interaction term in spin models) smoothes out the sharp spatial boundary 
between the phases predicted by the LHA. 

Even more interestingly, these similarities are not only qualitative . The distance from the spatial critical point 
where the switch between the two regimes takes place, x, can be determined by looking at the length scales relevant 
for the spatial quench. We compare the healing length 

(3) 



to the length scale on which the spatial critical point is approached: e/\de/dx\ assuming that e — \q — qd changes in 
x-direction only. This brings us to the spatial analog of H]) 

^° ' (4) 



Iq-Qcl" \de/dx\' 
To solve dH) we define — Ag^ and get 

The results © and (O are analogous. Indeed, apart from a slight difference in the scaling exponents - due to the 
absence of z in ([5]) - i maps onto x when time scales h/ Aq and tq are replaced by length scales and Aq, respectively. 
This shows striking parallels between quenches in time and space. It is also instructive to note that the same general 
result, Eq. can be obtained in a different way from the scaling theory (iTj . 



II. MODEL 



In the following, we will study a QPT in space in a ferromagnetic spin-1 Bosc-Einstein condensate. We consider 
untrapped clouds: atoms in a box. Such a system can be realized in an optical box trap [20| . Assuming that the 
condensate is placed in the magnetic field B{r) aligned in the z direction, its mean-field energy functional reads [2l| 

a 

where Fa=x,y,z is the spin-1 matrix. The strength of spin independent interactions is Cq = 47rfi,^(ao + 2a2)/3M > 0, 
while the strength of spin dependent interactions reads ci = 47r?i^(a2 — ao)/3M < 0. The constant as is the s-wave 
scattering length in the total spin S channel (ag = 101. Sa^, 02 = 100. 40^), and M is the atom mass. The prefactors 
of linear and quadratic Zeeman shifts are given by P = iibB{y)/2 and Q = pi^B'^{r)/AEhf, respectively. There E^f 
is the hyperfine splitting energy. 

The wave function has three condensate components, ipmj corresponding to m = 0, ±1 projections of spin-1 onto 
the magnetic field: ^'"^ = ("01, V'o, V'-i)- It is normalized as / dr^^xji — jy^ where N is the total number of atoms. 
The condensate magnetization reads 

/a = (*l^a|^'), a^x,y,z. 
It is convenient to define a dimensionless parameter 

g(r) =Q(r)/n|ci|, 

where n « N/V is the condensate density {V is the system volume). We consider below setups where n fluctuates 
negligibly in space. The critical point corresponds to qc — 2. 
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FIG. 1: A typical picture of the condensate magnetization in a spatial quench (the black line). The dashed red line shows the 
magnetization in the local homogeneous approximation, where the healing length is assumed to be negligibly small compared 
to the length scale of variations imposed on the condensate. The width of the crossover region is proportional to x to left and 
right of the spatial critical point at Xc = 2Aq [q{xc) = 2 ([7])]. This numerical simulation is for Xq/Xo = 3/2 - see Appendix iBl 
for numerical details. 



III. QUENCH IN A FERROMAGNETIC SPIN-1 BOSE-EINSTEIN CONDENSATE 

Similarly as in our work on time-dependent quench [ill, we assume that f dxfz — 0. Experimentally, such a 
constraint can be achieved b y p utting all atoms into an equal mixture of m = ±1 magnetic sublevels and letting the 
sample relax to equilibrium [2ll |. The spin conservation ensures that the final and initial states will have the same 
total fz magnetization: / drfz = 0. Additionally, we consider V'±i,o > in the ground state: a condition that can be 
always imposed because £[(|V'i|e'^\ |i/'o|e'^'', |V'-i 16*-^"^] > |i/'o|j IV'-il)]- This sets /y(r) = 0. Moreover, our 

numerical simulations indicate that /z(r) in the spatial quench considered here. Thus we investigate the /a;(r) 
condensate magnetization only. 

For simplicity, we assume that 

9(r) = ^, (7) 

and forgetting about the gradient term in © - i.e., using our LHA - the ground state phase diagram can be 
sketched. The part of the condensate exposed to < g(r) < 2 is in the ground state of the broken-symmetry 
phase. There V±i (r) = sfn^Xj^ — <7(r)/8 and i\j^^{v) = ^/nyJ\/2 + g(r)/4. The condensate is magnetized: fx(j) = 
ny/l — g(r)2/4, while /j,.2(r) = 0. This ground state breaks rotational symmetry on the (cc, y) plane present in energy 
functional ([6]). Parts of the condensate exposed to q{r) > 2 are in the polar phase ground state with ^^f{r) = and 
ip^^lr) — y/n, which implies fx,y,z{Y) — 0. The dependence of the condensate magnetization on x and q in the LHA 
is depicted with the red dashed line on Fig. [TJ There one condensate accommodates the two phases. 

The inclusion of the gradient term in the energy functional (|6]) prohibits singularities in the first derivative of the 
wave function. The smooth crossover region replaces the sharp phase boundary predicted by the LHA (Fig. [T|). Its 
size is proportional to x (O, as will be carefully discussed below. 

In the simplest setup q{v) dependence ([7]) is achieved by using the static, inhomogeneous, magnetic field to drive the 
system between the two phases. That would correspond to the following spatial variation of the Zeeman coefficients 
across a typical ^''Rb condensate: AQ '-^ n\ci\ and AP ~ 10^n|ci|. The latter result is obtained by using P = 
/i X 70 X 10^ X B HzG-\ Q = hx 70B^ HzC'^ 0, and n\ci\/h = 9.77 (evaluated at a peak density of the Berkeley 
experiment i). As can be found numerically, such an abrupt P- variation makes the LHA phase diagram qualitatively 
incorrect. 

To avoid these complications we propose to expose the condensate to the inhomogeneous magnetic field B pointing 
in the z-direction and changing in time faster than the characteristic time scales for the condensate dynamics. When 
{B{x,t))t — 0, where (• • •)( denotes time average, the linear Zeeman shift will be wiped out (which we assume from 
now on), while the quadratic term will be position-dependent only: Q{x) ~ {B{x,t)'^)t. Thus, we are still dealing 
with a phase transition in space. 

To make use of the general predictions worked out in ([5]), we have to determine the value of the critical exponent 
V and the prefactor from the generic expression for the healing length ([3]). We do it by comparing ([3]) to the 
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FIG. 2: Illustration of various scalings. The black dots: numerics showing the condensate magnetization at the spatial critical 
point = fx(xc = 2Aq). The green line presents the curve ~ ^q^^^- The red crosses and circles show numerics for x on 
the broken-symmetry side of the QPT. The blue line presents power law exp( — 1.7984)(Aq/Ao)'' '^°*'^, that was obtained from 
the fit Inx/Xo = -1.7984 ± 0.0005 + (0.3082 ± 0.0004) InAg/Ao, where Ao is a unit of length (see Appendix H]) . The fit was 
done to data in the range of Aq/Aq £ (1/2,8) - red crosses. Red circles show numerical data that was not used for this fit. 
The determination of x is illustrated on the inset showing the data for Xq/Xq = 3/2: x is the distance from the critical point 
where the relative departures of the condensate magnetization from the LHA result, f^^^ — ri^JX ~ q{x)'^/4:, start to exceed 
10% (other thresholds, 1% — 10%, give a similar result). Note logarithmic scale on both axis. 



expression for a divergent healing length near the critical point (jA3p . which we have derived in Appendix |X1 We get 
that = 1/2 on both sides of the QPT. Introducing a constant 

= ^hy2Mn\cil (8) 

we also see from (|A3I) that equals in the polar phase and ^s/ ■\/2 — 2|ci|/co in the broken-symmetry phase. 
Coming back to ^ we get that the size of the crossover region on either the polar or the broken-symmetry side is 

i^Co^^A^/^ (9) 

This expression can be verified by looking at experimentally measurable quantity: the ground state magnetization of 
the condensate @. 

We start our discussion on the polar side and linearize the wave-function around polar phase ground state: ■00 — 
\/n{l+S^o)7 ip±i = VnSilj±i, with Sip^ ^ 1- This leads to fx = V2n{6ipi+S^-i)+0{SilJm)- Using position-dependent 
q from ([7]) and linearized version of (|A1[) we arrive at 

4'VV.= (^-2)/., (10) 

which is solved by 

Ur) = ff'M (Ax/(AqC,2)1/3) /Ai(0). (11) 

There is the condensate magnetization at the spatial critical point (it depends on Xq only). Ax — x — 2Xq > is 
the distance from the spatial critical point, and Ai is an Airy function - the only nondivergent solution of (|10D . This 
solution rigorously shows that decay of the magnetization, /^^(r)//"'*, takes place on a length scale (Aq^^)^/^ in full 
agreement with the strikingly simple scaling result ([9]). Therefore, in the polar phase, the magnetization approaches 
the LHA value, fx = 0, at a distance x from the critical point. 

On the broken-symmetry side of the QPT, we refer to numerics (Appendix |B]) to verify accuracy of ([HI) . In the two 
extreme limits, very large and very small Xq the theory does not work well. Indeed, using ^ we see that the spatial 
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extent of the broken-symmetry phase is 2Aq (Fig. [T]). Thus, we expect that x <C 2Aq, i.e., Xq ^ to see well 

the crossover region on the broken-symmetry side. In the other limit, large Ag, we have to take into account that the 
system size in numerical simulations is finite, say I. There we need to have Ag ^ 1/2 so that there will be a spatial 
point in our system where q{x) — 2, and the finite size corrections coming from proximity of the system boundary to 
the spatial critical point will be negligible. These two estimates tell us that Ag in our system shall be much larger 
than O.lAo and much smaller than 39Ao (see Appendix [B] for numerical value of /, ~ ^s/V2 and the unit of length 
Aq relevant for our simulations). Increasing (decreasing) the lower (upper) limit on Ag by a factor of five we fit a 
power law to Ag/Ao G (1/2,8). We get from a fit that Inx/Ao = -1.7984 ± 0.0005 -I- (0.3082 ± 0.0004) In Ag/Ao - see 
Fig. [2l The numerical data for x clearly departs from the fitted line for Ag/Ap > 10. Departures for Ag/Ap ;$ 1/2 are 
much less pronounced. Therefore, we get approximately that x ^ ^^q'^ i which is in good qualitative agreement with 

1/3 

the predicted scaling, i.e., Ag . We expect that better agreement can be obtained for large systems where we can 

— 1/3 

explore larger Ag's for which crossover region is closer to the critical point {x/Xq ^ Xq ). More precisely, by taking 
the limit of N, /, Ag — > oo at = const (i.e., N/l = const) we expect that the scaling (|9]) will be fully recovered from 
numerical simulations. 

Now let's focus on the magnetization at the critical point. From (llip we know that it scales in the same way as the 
condensate magnetization at the border of the crossover region in the polar phase {x^ — 2Ag -I- x): 

Assuming that similar scaling relation will also hold between and fx at a;_ — 2Xq — x, i.e., at the broken-symmetry 
phase border of the crossover region where 

/,(x„)/n « v/l-9(^-)V4 ~ ^/s^ - Ag'/^ 

we get that /^"* - ^g^^^- As depicted in Fig. H we have indeed a robust ^ -^g^^^ scaling for Ag/Ao in the large 
range of 10~^ to 10. It is quite remarkable that the above intuitive prediction of the f^"^ scaling matches numerics 
so well. It suggests that magnetization at the spatial critical point might be a robust observable for studies of spatial 
quenches. Similar observation was made for condensate magnetization at the critical point during a temporal quench 
from the broken-symmetry to the polar phase p^ . 

All these results are in qualitative agreement with our expectations: in the limit of large Ag - very slow spatial 
driving - the condensate ground state approaches the LHA. Indeed, both the size of the crossover region in the 
q-parameter space, x/Xq, and the condensate magnetization at the critical point, /"'*, go to zero as Ag increases. 

The measurements of the dependence of the size of crossover region on Ag should allow for the first experimental 
determination of the scaling exponent v in a, ferromagnetic spin-1 Bose-Einstein condensate. Thus, our considerations 
lead to the new way of investigating the critical region of phase transitions. Another scheme that can be used to 
determine the exponent v was explored experimentally by Esslinger's group in the context of the classical phase 
transition from a normal gas to a Bose-Einstein condensate [2^ . The setup involved in this experiment is quite 
complicated as it requires presence of a high finesse optical cavity for detection purposes. We expect that the 
approach proposed by us should be easier to implement in the context of quantum phase transitions. Moreover, 
according to ([5]) our scheme shall be also ready for experimental exploration in the whole zoo of other physical 
systems undergoing a QPT. Finally, we would like to stress that applicability of Eq. ([5]) is not limited to mean-field 

1 /3 

theories: any experimental departures from the Ag scaling (apart from finite size corrections) should indicate that 
mean-field value of the critical exponent v does not describe the condensate properly. In such a case the correct value 
of the v exponent should be easily extracted from experimental data with the help of Eq. ([5]). 



IV. SUMMARY 

We have explored physics of the QPT in space in a ferromagnetic spin-1 Bose-Einstein condensate |6j: the most 
fiexible system for studies of QPTs in space and time. Our scaling results explain how singularities of the critical 
point affect the ground state magnetization of the condensate. They also suggest a new way for the measurement of 
the critical exponent v. These findings are generally applicable to all systems undergoing a second order QPT. The 
quest for the full understanding of a spatial quench opens up a prospect of interdisciplinary studies similarly as time 
quenches have done to date Q. 



6 



V. ACKNOWLEDGMENTS 

We acknowledge the support of the U.S. Department of Energy through the LANL/LDRD Program. 

APPENDIX A: HEALING LENGTH OF A FERROMAGNETIC SPIN-1 CONDENSATE 

The key ingredient of our theory is the divergent heahng length ([3]). Below we derive it from mean-field equations 

- ^V^V-ii +co(|?/'o|' + + \^-if)ip±i + cidi^oP + iV'iiP - lV'Til')V'±i +ciV'oV';i + W±i = AiV'ii, 

"I^V^Vo + codVoP + Wi? + IV^-iH^o + ci(|0ip + iV^-iH^o + 2ciVoV'iV'-i = mV-o, (A1) 

that come from minimization oi £ — ji J drj^m IV'mP without the linear Zeeman term ^ - see Sec. IIIII for the 
explanation of why we remove it from our considerations. 

The healing length is a typical length scale over which a local perturbation of the wave-function gets forgotten (here 
we will have three characteristic length scales, but we will identify the leading one relevant for a long distance healing 
process). To find it, we assume the constant parameter q — Q/n\ci\ across the condensate, and linearize mean-field 
equations (|Aip . We write the wave-function as tpm — "0™^ + \/nS^jjm with m — 0, ±1 and tp^ being the ground 
state solution. In the broken-symmetry phase we have tp±.f (r) — ^Jn^JXj^ — q{r)/8 and tpQ^ (r) — + q{r)/4:, 

while in the polar phase, '0±i(r) = and V'o^^(r) = -y/n. We use here the freedom to work with ■00, ±i > 0, which is 
explained in Sec. IIIII 

Linearized equations take the form 

^^V^J^™ = ^5„,fe<5Vfe, (A2) 



where = \/fi.^/2Mn|ci| and S stands for 

Co 



4|ci 

in the broken-symmetry phase, and for 



2 - qu_ (4 + 2(j)m_ 2 - qu+ 

{2-q)u^ 4 + 2q (2 - | , m± = 1 ± 2|ci|/co 

2 - qu+ (4 -f 2q)u- 2 - qu- 



q-\ -1 
2co/|ci| 
-1 q-\ 

in the polar phase. The factor u± above is very close to unity as |ci|/co ~ 1/216.1 ^ 1 for ^^Rb atoms considered 
here. 

To proceed we diagonalize matrix S by the transformation S = CDC~^, where D = diag(ri, Fa) is a diagonal 
matrix made of eigenvalues of matrix S. In the polar phase we find 

r„i = q - 2,g, 2co/|ci|, 

while in the broken-symmetry phase we get 

Tm = co/\ci\ T V(co/ci)2 - (4 - g2)(co/|ci| - l),q. 
After defining Sipm — X)n ^mli^i'n and some elementary algebra we get 



r- 7- l/L ~ T IIL ■ " T lit \ C I TP 

™ V -L ■ 

where the form of the function / depends on system dimensionality (ID, 2D, or 3D), while r is "attached" to the 
point in space where the perturbation is imposed on the wave-function. This solution implies that the eigenvalue 
vanishing at the critical point, Fmin, provides the leading contribution to the long-distance healing process near the 
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critical point (g « 2). Indeed, a simple calculation shows that the wave function will forget about the perturbation 
imposed on it {S-ipm, Sipm « 0) at a distance scaling as ^s/V^min- Thus, the divergent healing length equals £,s/V^min- 
Around the critical point we find it to be 



(^broken £^ t^polar ( 

" V2-2|ci|/cok-2|V2' « "|g-2|i/2 ^ ^ 

on the broken-symmetry (q < 2) and polar (q > 2) sides of the QPT, respectively. The expansion around q = 2 



was applied to obtain ^i^™kcn^ j have approximated the eigenvalue co/|ci| — \/{cq/ci)'^ — (4 — q'^){co/\ci \ — 1) 

by (2 — 2|ci|/co)(2 — q) near the critical point. These expressions fix the value of the mean-field critical exponent v 
for our model: ly — 1/2 (compare (jA3p to ([3])). We conclude by showing another simplification of the full expression 
for the divergent healing length on the broken symmetry side. Taking the limit of |ci|/co — s- 0, well justified for ®^Rb, 
we get 

^broken 



without restricting ourselves to q « 2. 



APPENDIX B: DETAILS OF NUMERICAL SIMULATIONS 

Our analytical predictions are valid in any dimensional system, but we test them numerically in a ID model. Our 
numerics is based on conjugate gradient minimization of the ID energy functional 

+ + l ^ {MFo.\^r+qfi\c^mF!\i'), (Bi) 

a—x^y^z 

obtained through dimensional reduction of the full 3D energy functional ^ without the linear Zeeman term - see Sec. 
mil for the explanation of why we remove it from our considerations. Moreover, all the components of the above energy 
functional are dimensionless, and chosen to match qualitatively the experiment of the Berkeley group [6J. Derivation 
of (|B1|) is explained in detail in [ll|: please account for different normalization of the wave- function between this 
paper and jllj. 

In short, we have N, number of atoms, equal to 2 x 10^, Nci = —629, and cq — 216.1|ci|. The unit of length 
used through the paper, Aq, equals b.lfim. The system (box) size is I = l/Xo = 78, i.e., about OAmm {I x Aq). The 
coordinate x runs from to I. The density, h, equals N/l, while J dx^^^ — N. The "prefactor" of the healing 

length, ([HI), in these dimensionless quantities is yl/2N\ci\ « 1/4, i.e., about 1.3/im (Ao/4). 
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